## Explaining Attitudes Toward Immigration: The Role of Regional Context and Individual Predispositions
## R script to reproduce analyses in the article & online appendix
## Johannes Karreth
## Last updated 2020/08/21

## NOTE: I updated this script on 2020/08/21 after noticing that the update
## to R 4.0 and the changed default for stringsAsFactors breaks the original 
## code.

## Please see "readme.txt" for a list of all data files that need to be 
## put in your working directory

####################################
## Set working directory
####################################

# setwd("...")

####################################
## Load packages
####################################

library(foreign)
library(reshape2)
library(lme4)
library(effects)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(texreg)
library(polycor)
library(maptools)
library(dplyr)

####################################
## My R session info before continuing this file
####################################

# R version 4.0.2 (2020-06-22)
# Platform: x86_64-apple-darwin17.0 (64-bit)
# Running under: macOS Catalina 10.15.6
# 
# Matrix products: default
# BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
# 
# locale:
#   [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
#   [1] dplyr_1.0.1    maptools_1.0-1 sp_1.4-2       polycor_0.7-10 texreg_1.37.5  gridExtra_2.3 
# [7] ggthemes_4.2.0 ggplot2_3.3.2  effects_4.2-0  carData_3.0-4  lme4_1.1-23    Matrix_1.2-18 
# [13] reshape2_1.4.4 foreign_0.8-80
# 
# loaded via a namespace (and not attached):
#   [1] statmod_1.4.34   tidyselect_1.1.0 xfun_0.16        purrr_0.3.4      mitools_2.4     
# [6] splines_4.0.2    lattice_0.20-41  colorspace_1.4-1 vctrs_0.3.2      generics_0.0.2  
# [11] htmltools_0.5.0  yaml_2.2.1       survival_3.2-3   rlang_0.4.7      nloptr_1.2.2.2  
# [16] pillar_1.4.6     glue_1.4.1       withr_2.2.0      DBI_1.1.0        lifecycle_0.2.0 
# [21] plyr_1.8.6       stringr_1.4.0    munsell_0.5.0    gtable_0.3.0     evaluate_0.14   
# [26] knitr_1.29       Rcpp_1.0.5.2     scales_1.1.1     digest_0.6.25    stringi_1.4.6   
# [31] insight_0.9.0    survey_4.0       grid_4.0.2       tools_4.0.2      magrittr_1.5    
# [36] tibble_3.0.3     crayon_1.3.4     pkgconfig_2.0.3  MASS_7.3-51.6    ellipsis_0.3.1  
# [41] httr_1.4.2       minqa_1.2.4      rmarkdown_2.3    rstudioapi_0.11  R6_2.4.1        
# [46] boot_1.3-25      nnet_7.3-14      nlme_3.1-148     compiler_4.0.2   

####################################
## Figures 1 & 3: Maps
####################################

shape <- rgdal::readOGR(dsn = "de_at_ch")
shape@data$id <- shape@data$ccbu
shape.points <- fortify(shape, region = "ccbu")
shape.df <- left_join(x = shape.points, y = shape@data, by = "id")

map.dat <- read.csv("mapdata.csv", stringsAsFactors = TRUE)
map.dat$id <- map.dat$ccbu
regions.mapdata <- left_join(x = shape.df, y = map.dat, by = "id")

p.ESSimdfetn2 <- ggplot(regions.mapdata) + aes(long, lat, group = group, fill = cut(ESSimdfetn2, breaks = 8, dig.lab = 0)) + geom_polygon() + coord_equal() + geom_path(color = "black") + scale_fill_brewer(palette = "Greys", type = "seq") + scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) + xlab("") + ylab("") + theme_tufte(base_family = "sans", base_size = 16) + theme(legend.position=c(0.95,1), legend.justification=c(1,1)) + labs(fill = "") + ggtitle("Allow many immigrants of \ndifferent race/ethnic group?")

p.ESSimpcntr2 <- ggplot(regions.mapdata) + aes(long, lat, group = group, fill = cut(ESSimpcntr2, breaks = 8, dig.lab = 0)) + geom_polygon() + coord_equal() + geom_path(color = "black") + scale_fill_brewer(palette = "Greys", type = "seq") + scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) + xlab("") + ylab("") + theme_tufte(base_family = "sans", base_size = 16) + theme(legend.position=c(0.95,1), legend.justification=c(1,1)) + labs(fill = "") + ggtitle("Allow many immigrants from \npoorer countries outside Europe?")

# Print figures

pdf(file = "attitudes_map.pdf", width = 20, height = 10)
grid.arrange(p.ESSimdfetn2, p.ESSimpcntr2, ncol = 2)
dev.off()

p.fnrbch91 <- ggplot(regions.mapdata) + aes(long, lat, group = group, fill = cut(fnrbch91, breaks = 9, dig.lab = 0)) + geom_polygon() + coord_equal() + geom_path(color = "black") + scale_fill_brewer(palette = "Greys", type = "seq") + scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) + xlab("") + ylab("") + theme_tufte(base_family = "sans", base_size = 16) + theme(legend.position=c(0.95,1), legend.justification=c(1,1)) + labs(fill = "") + ggtitle(expression("Change in % Foreign \npopulation (1991-2008)"))

p.nwepr2008 <- ggplot(regions.mapdata) + aes(long, lat, group = group, fill = cut(nwepr2008, breaks = 11, dig.lab = 0)) + geom_polygon() + coord_equal() + geom_path(color = "black") + scale_fill_brewer(palette = "Greys", type = "seq") + scale_x_continuous(breaks = NULL) + scale_y_continuous(breaks = NULL) + xlab("") + ylab("") + theme_tufte(base_family = "sans", base_size = 16) + theme(legend.position=c(0.95,1), legend.justification=c(1,1)) + labs(fill = "") + ggtitle("% Non-Western immigrants (2008)")

pdf(file = "migrants_map.pdf", width = 20, height = 10)
grid.arrange(p.fnrbch91, p.nwepr2008, ncol = 2)
dev.off()

# For conversion to EPS, use the pdf2ps and then ps2eps programs via the command line
# Example
# pdf2ps "migrants_map.pdf"
# ps2eps "migrants_map.ps"

####################################
## Figure 2: Time trends of % foreign population
####################################

foreigntrend.dat <- read.csv("regiondata.csv", stringsAsFactors = TRUE)
names(foreigntrend.dat)

foreigntrend.dat$fnrb_i <- foreigntrend.dat$fnrb_i*100
foreigntrend.dat$fnrb <- foreigntrend.dat$fnrb*100
foreigntrend.dat$uniqid <- as.factor <- foreigntrend.dat$uniqid

foreigntrend.dat <- foreigntrend.dat[foreigntrend.dat$year > 1990, ]

ts <- ggplot(data = foreigntrend.dat, mapping = aes(x = year, y = fnrb_i, group = uniqid)) + geom_point() + geom_line() + geom_smooth(mapping = aes(x = year, y = fnrb, group = NULL)) + facet_wrap(~ country, scales = "free") + xlab("Year") + ylab("% Foreign Population") + theme_bw()

pdf("fnr_ts.pdf", width = 8, height = 4)
print(ts)
dev.off()

####################################
## Country-level analyses
####################################

####################################
## Prep data
####################################

country.dat <- read.csv("ESS1_country.csv", stringsAsFactors = TRUE)
country.dat <- country.dat[country.dat$ESStnat == 1, ]

country.dat$ESSagegroup2 <- as.numeric(country.dat$ESSagegroup)

table(country.dat$ESSimdfetn2)
table(country.dat$ESSimpcntr2)

country.dat$ESSimdfetn2.std <- scale(x = country.dat$ESSimdfetn2, center = TRUE, scale = TRUE)[, 1]
country.dat$ESSimpcntr2.std <- scale(x = country.dat$ESSimpcntr2, center = TRUE, scale = TRUE)[, 1]

## Models C1 - C8

c.mod1 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat)
c.mod2 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale*fnr_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + (ESSlrscale | ccode), data = country.dat)
c.mod3 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat)
c.mod4 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + (ESSlrscale | ccode), data = country.dat)

c.mod5 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat)
c.mod6 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale*fnr_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + (ESSlrscale | ccode), data = country.dat)
c.mod7 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat)
c.mod8 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup2 + ESSlrscale*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + (ESSlrscale | ccode), data = country.dat)

screenreg(list(c.mod1, c.mod2, c.mod3, c.mod4, c.mod5, c.mod6, c.mod7, c.mod8), digits = 3)

ctab1 <- texreg(list(c.mod1, c.mod2, c.mod3, c.mod4), 
               file = "cmods1.tex", 
               stars = 0.05, 
               custom.model.names = c("C1", "C2", "C3", "C4"), 
               custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
               custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
               digits = 3, 
               leading.zero = TRUE, 
               omit.coef = "(Intercept)", 
               center = TRUE, 
               caption = "Country-level results. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
               caption.above = TRUE, 
               label = "tab:countryregs1", 
               dcolumn = TRUE, 
               booktabs = TRUE, 
               use.packages = FALSE, 
               table = TRUE)

ctab1 <- htmlreg(list(c.mod1, c.mod2, c.mod3, c.mod4), 
                file = "cmods2.doc", 
                stars = 0.05, 
                custom.model.names = c("C1", "C2", "C3", "C4"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Country-level results. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE
                )

ctab2 <- texreg(list(c.mod5, c.mod6, c.mod7, c.mod8), 
                file = "cmods2.tex", 
                stars = 0.05, 
                custom.model.names = c("C5", "C6", "C7", "C8"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Country-level results. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:countryregs2", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

ctab2 <- htmlreg(list(c.mod5, c.mod6, c.mod7, c.mod8), 
                file = "cmods2.doc", 
                stars = 0.05, 
                custom.model.names = c("C5", "C6", "C7", "C8"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Country-level results. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE)

####################################
## Figure 4: Conditional effects (country-level)
####################################

c.eff2 <- effect("ESSlrscale*fnr_ch91", c.mod2, xlevels = list(ESSlrscale = c(1, 6, 11), fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
c.eff4 <- effect("ESSlrscale*nwepr", c.mod4, xlevels = list(ESSlrscale = c(1, 6, 11), nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)
c.eff6 <- effect("ESSlrscale*fnr_ch91", c.mod6, xlevels = list(ESSlrscale = c(1, 6, 11), fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90)
c.eff8 <- effect("ESSlrscale*nwepr", c.mod8, xlevels = list(ESSlrscale = c(1, 6, 11), nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)

c.eff2.dat <- cbind(c.eff2$fit, c.eff2$x, c.eff2$lower, c.eff2$upper); names(c.eff2.dat) <- c("fit", "lr", "fnr_ch91", "lower", "upper")
c.eff2.dat$lr2 <- as.factor(c.eff2.dat$lr); levels(c.eff2.dat$lr2) <- c("Left", "Center", "Right")
c.eff4.dat <- cbind(c.eff4$fit, c.eff4$x, c.eff4$lower, c.eff4$upper); names(c.eff4.dat) <- c("fit", "lr", "nwepr", "lower", "upper")
c.eff4.dat$lr2 <- as.factor(c.eff4.dat$lr); levels(c.eff4.dat$lr2) <- c("Left", "Center", "Right")
c.eff6.dat <- cbind(c.eff6$fit, c.eff6$x, c.eff6$lower, c.eff6$upper); names(c.eff6.dat) <- c("fit", "lr", "fnr_ch91", "lower", "upper")
c.eff6.dat$lr2 <- as.factor(c.eff6.dat$lr); levels(c.eff6.dat$lr2) <- c("Left", "Center", "Right")
c.eff8.dat <- cbind(c.eff8$fit, c.eff8$x, c.eff8$lower, c.eff8$upper); names(c.eff8.dat) <- c("fit", "lr", "nwepr", "lower", "upper")
c.eff8.dat$lr2 <- as.factor(c.eff8.dat$lr); levels(c.eff8.dat$lr2) <- c("Left", "Center", "Right")

rug.dat <- data.frame(fnr_ch91 = country.dat$fnr_ch91, nwepr = country.dat$nwepr)

c.eff2.p <- ggplot(dat = c.eff2.dat, aes(x = fnr_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnr_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (C2)") + geom_rug(data = rug.dat, aes(x = jitter(fnr_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

c.eff4.p <- ggplot(dat = c.eff4.dat, aes(x = nwepr, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (C4)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

c.eff6.p <- ggplot(dat = c.eff6.dat, aes(x = fnr_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnr_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (C6)") + geom_rug(data = rug.dat, aes(x = jitter(fnr_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue")) # + theme(legend.position = c(0.1, 0.2))

c.eff8.p <- ggplot(dat = c.eff8.dat, aes(x = nwepr, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (C8)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

pdf("int_country.pdf", width = 8, height = 10)
grid.arrange(c.eff2.p, c.eff4.p, c.eff6.p, c.eff8.p, ncol = 2)
dev.off()

####################################
## Split samples L/C/R (country-level): C1-Left - C7-Right
####################################

c.mod1.left.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale < 6, ])
c.mod1.center.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale == 6, ])
c.mod1.right.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale > 6, ])
c.mod3.left.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale < 6, ])
c.mod3.center.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale == 6, ])
c.mod3.right.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale > 6, ])

c.mod5.left.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale < 6, ])
c.mod5.center.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale == 6, ])
c.mod5.right.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + fnr_ch91 + (1 | ccode), data = country.dat[country.dat$ESSlrscale > 6, ])
c.mod7.left.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale < 6, ])
c.mod7.center.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale == 6, ])
c.mod7.right.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup2 + ESSedulvla2 + ESShinc2 + unemp + fnr1991 + nwepr + (1 | ccode), data = country.dat[country.dat$ESSlrscale > 6, ])


screenreg(list(c.mod1.left.tab, c.mod1.center.tab, c.mod1.right.tab, c.mod3.left.tab, c.mod3.center.tab, c.mod3.right.tab), digits = 3)

screenreg(list(c.mod5.left.tab, c.mod5.center.tab, c.mod5.right.tab, c.mod7.left.tab, c.mod7.center.tab, c.mod7.right.tab), digits = 3)

# Tables

ctab1.lcr <- texreg(list(c.mod1.left.tab, c.mod1.center.tab, c.mod1.right.tab, c.mod3.left.tab, c.mod3.center.tab, c.mod3.right.tab), 
                    file = "countryregs1_lcr.tex", 
                    stars = 0.05, 
                    custom.model.names = c("C1-Left", "C1-Center", "C1-Right", "C3-Left", "C3-Center", "C3-Right"), 
                    custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "% Non-Western Foreigners"),
                    custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
                    digits = 3, 
                    leading.zero = TRUE, 
                    omit.coef = "(Intercept)", 
                    center = TRUE, 
                    caption = "Country-level results, samples split in respondents on the political Left, Center, and Right. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                    caption.above = TRUE, 
                    label = "tab:countryregs1_lcr", 
                    dcolumn = TRUE, 
                    booktabs = TRUE, 
                    use.packages = FALSE, 
                    table = TRUE)

ctab2.lcr <- texreg(list(c.mod5.left.tab, c.mod5.center.tab, c.mod5.right.tab, c.mod7.left.tab, c.mod7.center.tab, c.mod7.right.tab), 
                    file = "countryregs2_lcr.tex", 
                    stars = 0.05, 
                    custom.model.names = c("C5-Left", "C5-Center", "C5-Right", "C7-Left", "C7-Center", "C7-Right"), 
                    custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "% Non-Western Foreigners"),
                    custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for countries (standard errors in parentheses).", 
                    digits = 3, 
                    leading.zero = TRUE, 
                    omit.coef = "(Intercept)", 
                    center = TRUE, 
                    caption = "Country-level results, samples split in respondents on the political Left, Center, and Right. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                    caption.above = TRUE, 
                    label = "tab:countryregs2_lcr", 
                    dcolumn = TRUE, 
                    booktabs = TRUE, 
                    use.packages = FALSE, 
                    table = TRUE)

## Plots

c.mod1.left.eff <- effect("fnr_ch91", c.mod1.left.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
c.mod1.center.eff <- effect("fnr_ch91", c.mod1.center.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90)
c.mod1.right.eff <- effect("fnr_ch91", c.mod1.right.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90)

c.mod3.left.eff <- effect("nwepr", c.mod3.left.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
c.mod3.center.eff <- effect("nwepr", c.mod3.center.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)
c.mod3.right.eff <- effect("nwepr", c.mod3.right.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)

c.mod5.left.eff <- effect("fnr_ch91", c.mod5.left.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
c.mod5.center.eff <- effect("fnr_ch91", c.mod5.center.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90)
c.mod5.right.eff <- effect("fnr_ch91", c.mod5.right.tab, xlevels = list(fnr_ch91 = sort(unique(country.dat$fnr_ch91))), confidence.level = 0.90)

c.mod7.left.eff <- effect("nwepr", c.mod7.left.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
c.mod7.center.eff <- effect("nwepr", c.mod7.center.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)
c.mod7.right.eff <- effect("nwepr", c.mod7.right.tab, xlevels = list(nwepr = sort(unique(country.dat$nwepr))), confidence.level = 0.90)

c.mod1.left.eff.dat <- data.frame(c.mod1.left.eff$fit, c.mod1.left.eff$x, c.mod1.left.eff$lower, c.mod1.left.eff$upper, LR = "Left"); names(c.mod1.left.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod1.center.eff.dat <- data.frame(c.mod1.center.eff$fit, c.mod1.center.eff$x, c.mod1.center.eff$lower, c.mod1.center.eff$upper, LR = "Center"); names(c.mod1.center.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod1.right.eff.dat <- data.frame(c.mod1.right.eff$fit, c.mod1.right.eff$x, c.mod1.right.eff$lower, c.mod1.right.eff$upper, LR = "Right"); names(c.mod1.right.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod1.lcr.eff.dat <- rbind(c.mod1.left.eff.dat, c.mod1.center.eff.dat, c.mod1.right.eff.dat)

c.mod3.left.eff.dat <- data.frame(c.mod3.left.eff$fit, c.mod3.left.eff$x, c.mod3.left.eff$lower, c.mod3.left.eff$upper, LR = "Left"); names(c.mod3.left.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod3.center.eff.dat <- data.frame(c.mod3.center.eff$fit, c.mod3.center.eff$x, c.mod3.center.eff$lower, c.mod3.center.eff$upper, LR = "Center"); names(c.mod3.center.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod3.right.eff.dat <- data.frame(c.mod3.right.eff$fit, c.mod3.right.eff$x, c.mod3.right.eff$lower, c.mod3.right.eff$upper, LR = "Right"); names(c.mod3.right.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod3.lcr.eff.dat <- rbind(c.mod3.left.eff.dat, c.mod3.center.eff.dat, c.mod3.right.eff.dat)

c.mod5.left.eff.dat <- data.frame(c.mod5.left.eff$fit, c.mod5.left.eff$x, c.mod5.left.eff$lower, c.mod5.left.eff$upper, LR = "Left"); names(c.mod5.left.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod5.center.eff.dat <- data.frame(c.mod5.center.eff$fit, c.mod5.center.eff$x, c.mod5.center.eff$lower, c.mod5.center.eff$upper, LR = "Center"); names(c.mod5.center.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod5.right.eff.dat <- data.frame(c.mod5.right.eff$fit, c.mod5.right.eff$x, c.mod5.right.eff$lower, c.mod5.right.eff$upper, LR = "Right"); names(c.mod5.right.eff.dat) <- c("fit", "fnr_ch91", "lower", "upper", "LR")
c.mod5.lcr.eff.dat <- rbind(c.mod5.left.eff.dat, c.mod5.center.eff.dat, c.mod5.right.eff.dat)

c.mod7.left.eff.dat <- data.frame(c.mod7.left.eff$fit, c.mod7.left.eff$x, c.mod7.left.eff$lower, c.mod7.left.eff$upper, LR = "Left"); names(c.mod7.left.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod7.center.eff.dat <- data.frame(c.mod7.center.eff$fit, c.mod7.center.eff$x, c.mod7.center.eff$lower, c.mod7.center.eff$upper, LR = "Center"); names(c.mod7.center.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod7.right.eff.dat <- data.frame(c.mod7.right.eff$fit, c.mod7.right.eff$x, c.mod7.right.eff$lower, c.mod7.right.eff$upper, LR = "Right"); names(c.mod7.right.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
c.mod7.lcr.eff.dat <- rbind(c.mod7.left.eff.dat, c.mod7.center.eff.dat, c.mod7.right.eff.dat)

rug.dat <- data.frame(fnr_ch91 = country.dat$fnr_ch91, nwepr = country.dat$nwepr)

c.mod1.lcr.eff.p <- ggplot(dat = c.mod1.lcr.eff.dat, aes(x = fnr_ch91, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = fnr_ch91, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (C1-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(fnr_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

c.mod3.lcr.eff.p <- ggplot(dat = c.mod3.lcr.eff.dat, aes(x = nwepr, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("% Non−Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (C3-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

c.mod5.lcr.eff.p <- ggplot(dat = c.mod5.lcr.eff.dat, aes(x = fnr_ch91, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = fnr_ch91, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (C5-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(fnr_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

c.mod7.lcr.eff.p <- ggplot(dat = c.mod7.lcr.eff.dat, aes(x = nwepr, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("% Non−Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (C7-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

pdf("int_country_lcr.pdf", width = 12, height = 10)
grid.arrange(c.mod1.lcr.eff.p, c.mod3.lcr.eff.p, c.mod5.lcr.eff.p, c.mod7.lcr.eff.p, ncol = 2)
dev.off()



####################################
## Region-level analyses
####################################

####################################
## Prep data
####################################


region.dat <- read.csv("ESS1-4_de-at-ch_region.csv", stringsAsFactors = TRUE)
region.dat <- region.dat[region.dat$ESStnat == 1, ]

region.dat$ESSagegroup <- as.numeric(region.dat$ESSagegroup)
region.dat$uniqid2 <- as.numeric(region.dat$uniqid)

table(region.dat$ESSimdfetn2)
table(region.dat$ESSimpcntr2)

region.dat$ESSimdfetn2.std <- scale(x = region.dat$ESSimdfetn2, center = TRUE, scale = TRUE)[, 1]
region.dat$ESSimpcntr2.std <- scale(x = region.dat$ESSimpcntr2, center = TRUE, scale = TRUE)[, 1]
region.dat$ESSagegroup.std <- region.dat$ESSagegroup - 6
region.dat$ESSedulvla2.std <- region.dat$ESSedulvla2 - 2
region.dat$ESShinc2.std <- region.dat$ESShinc2 - 5

region.dat$year <- ifelse(region.dat$yr2004 == 1, 2004, NA)
region.dat$year <- ifelse(region.dat$yr2006 == 1, 2006, region.dat$year)
region.dat$year <- ifelse(region.dat$yr2008 == 1, 2008, region.dat$year)
region.dat$year <- ifelse(region.dat$yr2004 == 0 & region.dat$yr2006 == 0 & region.dat$yr2008 == 0, 2002, region.dat$year)
region.dat$uniqidyr <- paste(region.dat$uniqid, region.dat$year)
region.dat$uniqid2yr <- region.dat$uniqidyr

####################################
## Main models: R1 - R8
####################################

r.mod1.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod2.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod3.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod4.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod5.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod6.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod7.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod8.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)

screenreg(list(r.mod1.tab, r.mod2.tab, r.mod3.tab, r.mod4.tab, r.mod5.tab, r.mod6.tab, r.mod7.tab, r.mod8.tab), digits = 3)

rtab1 <- texreg(list(r.mod1.tab, r.mod2.tab, r.mod3.tab, r.mod4.tab), 
                file = "regionregs1.tex", 
                stars = 0.05, 
                custom.model.names = c("R1", "R2", "R3", "R4"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs1", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab1 <- htmlreg(list(r.mod1.tab, r.mod2.tab, r.mod3.tab, r.mod4.tab), 
                file = "regionregs1.doc", 
                stars = 0.05, 
                custom.model.names = c("R1", "R2", "R3", "R4"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE)

rtab2 <- texreg(list(r.mod5.tab, r.mod6.tab, r.mod7.tab, r.mod8.tab), 
                file = "/Users/johanneskarreth/Documents/Uni/1 - Papers/26 - Diversity by State/Data/Tables/regionregs2.tex", 
                stars = 0.05, 
                custom.model.names = c("R5", "R6", "R7", "R8"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs2", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab2 <- htmlreg(list(r.mod5.tab, r.mod6.tab, r.mod7.tab, r.mod8.tab), 
                file = "regionregs2.doc", 
                stars = 0.05, 
                custom.model.names = c("R5", "R6", "R7", "R8"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE)

# Does the interaction improve the models?
anova(r.mod1, r.mod2)
anova(r.mod3, r.mod4)
anova(r.mod5, r.mod6)
anova(r.mod7, r.mod8)

####################################
## Random effects for region-years: R1-RY - R8-RY
####################################

r.mod1.ryre <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2yr), data = region.dat)
r.mod2.ryre <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2yr), data = region.dat)
r.mod3.ryre <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2yr), data = region.dat)
r.mod4.ryre <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2yr), data = region.dat)
r.mod5.ryre <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2yr), data = region.dat)
r.mod6.ryre <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2yr), data = region.dat)
r.mod7.ryre <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2yr), data = region.dat)
r.mod8.ryre <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2yr), data = region.dat)

screenreg(list(r.mod1.ryre, r.mod2.ryre, r.mod3.ryre, r.mod4.ryre, r.mod5.ryre, r.mod6.ryre, r.mod7.ryre, r.mod8.ryre), digits = 3)

rtab1.ryre <- texreg(list(r.mod1.ryre, r.mod2.ryre, r.mod3.ryre, r.mod4.ryre), 
                file = "regionregs1ryre.tex", 
                stars = 0.05, 
                custom.model.names = c("R1-RY", "R2-RY", "R3-RY", "R4-RY"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for region-years (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs1ryre", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab2.ryre <- texreg(list(r.mod5.ryre, r.mod6.ryre, r.mod7.ryre, r.mod8.ryre), 
                file = "regionregs2ryre.tex", 
                stars = 0.05, 
                custom.model.names = c("R5-RY", "R6-RY", "R7-RY", "R8-RY"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for region-years (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs2ryre", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

####################################
## Fixed effects for survey years & countries: R1-FE - R8-FE
####################################

r.femod1.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91+ yr2004 + yr2006 + yr2008 + austria + switzerland + (1 | uniqid2), data = region.dat)
r.femod2.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + yr2004 + yr2006 + yr2008 + austria + switzerland + (ESSlrscale2 | uniqid2), data = region.dat)
r.femod3.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + yr2004 + yr2006 + yr2008 + austria + switzerland + (1 | uniqid2), data = region.dat)
r.femod4.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + yr2004 + yr2006 + yr2008 + austria + switzerland + (ESSlrscale2 | uniqid2), data = region.dat)
r.femod5.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91+ yr2004 + yr2006 + yr2008 + austria + switzerland + (1 | uniqid2), data = region.dat)
r.femod6.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + yr2004 + yr2006 + yr2008 + austria + switzerland + (ESSlrscale2 | uniqid2), data = region.dat)
r.femod7.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + yr2004 + yr2006 + yr2008 + austria + switzerland + (1 | uniqid2), data = region.dat)
r.femod8.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + yr2004 + yr2006 + yr2008 + austria + switzerland + (ESSlrscale2 | uniqid2), data = region.dat)

screenreg(list(r.femod1.tab, r.femod2.tab, r.femod3.tab, r.femod4.tab, r.femod5.tab, r.femod6.tab, r.femod7.tab, r.femod8.tab), digits = 3)

rfetab1 <- texreg(list(r.femod1.tab, r.femod2.tab, r.femod3.tab, r.femod4.tab), 
                file = "regionferegs1.tex", 
                stars = 0.05, 
                custom.model.names = c("R1-FE", "R2-FE", "R3-FE", "R4-FE"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Survey year: 2004", "Survey year: 2006", "Survey year: 2008", "Austria", "Switzerland", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results, including fixed effects for survey years and countries. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs1", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rfetab2 <- texreg(list(r.femod5.tab, r.femod6.tab, r.femod7.tab, r.femod8.tab), 
                file = "regionferegs2.tex", 
                stars = 0.05, 
                custom.model.names = c("R5-FE", "R6-FE", "R7-FE", "R8-FE"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Survey year: 2004", "Survey year: 2006", "Survey year: 2008", "Austria", "Switzerland", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results, including fixed effects for survey years and countries.. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs2", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

####################################
## Recode attitudes to binary: R1-B - R8-B
####################################

region.dat$ESSimdfetn2.bin <- with(region.dat, ifelse(ESSimdfetn2 < 3, 0, 1))
region.dat$ESSimpcntr2.bin <- with(region.dat, ifelse(ESSimpcntr2 < 3, 0, 1))

r.mod1.bin.tab <- glmer(ESSimdfetn2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat, family = "binomial")
r.mod2.bin.tab <- glmer(ESSimdfetn2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat, family = "binomial")
r.mod3.bin.tab <- glmer(ESSimdfetn2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat, family = "binomial")
r.mod4.bin.tab <- glmer(ESSimdfetn2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat, family = "binomial")
r.mod5.bin.tab <- glmer(ESSimpcntr2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat, family = "binomial")
r.mod6.bin.tab <- glmer(ESSimpcntr2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat, family = "binomial")
r.mod7.bin.tab <- glmer(ESSimpcntr2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat, family = "binomial")
r.mod8.bin.tab <- glmer(ESSimpcntr2.bin ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat, family = "binomial")

screenreg(list(r.mod1.bin.tab, r.mod2.bin.tab, r.mod3.bin.tab, r.mod4.bin.tab, r.mod5.bin.tab, r.mod6.bin.tab, r.mod7.bin.tab, r.mod8.bin.tab), digits = 3)

rtab1bin <- texreg(list(r.mod1.bin.tab, r.mod2.bin.tab, r.mod3.bin.tab, r.mod4.bin.tab), 
                file = "regionregs1bin.tex", 
                stars = 0.05, 
                custom.model.names = c("R1", "R2", "R3", "R4"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical logit models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Binary outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Positive answers recoded as 1, negative as 0.", 
                caption.above = TRUE, 
                label = "tab:regionregs1bin", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab2bin <- texreg(list(r.mod5.bin.tab, r.mod6.bin.tab, r.mod7.bin.tab, r.mod8.bin.tab), 
                file = "regionregs2bin.tex", 
                stars = 0.05, 
                custom.model.names = c("R5", "R6", "R7", "R8"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical logit models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results. Binary outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Positive responses recoded as 1, negative as 0.", 
                caption.above = TRUE, 
                label = "tab:regionregs2bin", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

## Figures: conditional effects

rbin.eff2 <- effect("ESSlrscale2*fnrb_ch91", r.mod2.bin.tab, xlevels = list(ESSlrscale2 = c(0, 5, 10), fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90); # plot(r.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
rbin.eff4 <- effect("ESSlrscale2*nwepr", r.mod4.bin.tab, xlevels = list(ESSlrscale2 = c(0, 5, 10), nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)
rbin.eff6 <- effect("ESSlrscale2*fnrb_ch91", r.mod6.bin.tab, xlevels = list(ESSlrscale2 = c(0, 5, 10), fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)
rbin.eff8 <- effect("ESSlrscale2*nwepr", r.mod8.bin.tab, xlevels = list(ESSlrscale2 = c(0, 5, 10), nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)

rbin.eff2.dat <- cbind(exp(rbin.eff2$fit) / (1 + exp(rbin.eff2$fit)), rbin.eff2$x, exp(rbin.eff2$lower) / (1 + exp(rbin.eff2$lower)), exp(rbin.eff2$upper) / (1 + exp(rbin.eff2$upper))); names(rbin.eff2.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
rbin.eff2.dat$lr2 <- as.factor(rbin.eff2.dat$lr); levels(rbin.eff2.dat$lr2) <- c("Left", "Center", "Right")
rbin.eff4.dat <- cbind(exp(rbin.eff4$fit) / (1 + exp(rbin.eff4$fit)), rbin.eff4$x, exp(rbin.eff4$lower) / (1 + exp(rbin.eff4$lower)), exp(rbin.eff4$upper) / (1 + exp(rbin.eff4$upper))); names(rbin.eff4.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
rbin.eff4.dat$lr2 <- as.factor(rbin.eff4.dat$lr); levels(rbin.eff4.dat$lr2) <- c("Left", "Center", "Right")
rbin.eff6.dat <- cbind(exp(rbin.eff6$fit) / (1 + exp(rbin.eff6$fit)), rbin.eff6$x, exp(rbin.eff6$lower) / (1 + exp(rbin.eff6$lower)), exp(rbin.eff6$upper) / (1 + exp(rbin.eff6$upper))); names(rbin.eff6.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
rbin.eff6.dat$lr2 <- as.factor(rbin.eff6.dat$lr); levels(rbin.eff6.dat$lr2) <- c("Left", "Center", "Right")
rbin.eff8.dat <- cbind(exp(rbin.eff8$fit) / (1 + exp(rbin.eff8$fit)), rbin.eff8$x, exp(rbin.eff8$lower) / (1 + exp(rbin.eff8$lower)), exp(rbin.eff8$upper) / (1 + exp(rbin.eff8$upper))); names(rbin.eff8.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
rbin.eff8.dat$lr2 <- as.factor(rbin.eff8.dat$lr); levels(rbin.eff8.dat$lr2) <- c("Left", "Center", "Right")

rug.dat <- data.frame(fnrb_ch91 = region.dat$fnrb_ch91, nwepr = region.dat$nwepr)

rbin.eff2.p <- ggplot(dat = rbin.eff2.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Pr(Positive acceptance)") + ggtitle("Ethnic acceptance (R2)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = 0.45, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = 0.15, label = "Right", colour = "blue")+ scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

rbin.eff4.p <- ggplot(dat = rbin.eff4.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Pr(Positive acceptance)") + ggtitle("Ethnic acceptance (R4)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = 0.45, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = 0.15, label = "Right", colour = "blue")+ scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

rbin.eff6.p <- ggplot(dat = rbin.eff6.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Pr(Positive acceptance)") + ggtitle("Economic acceptance (R6)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = 0.45, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = 0.15, label = "Right", colour = "blue")+ scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

rbin.eff8.p <- ggplot(dat = rbin.eff8.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Pr(Positive acceptance)") + ggtitle("Economic acceptance (R8)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = 0.45, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = 0.15, label = "Right", colour = "blue")+ scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

pdf("intbin_region.pdf", width = 8, height = 10)
grid.arrange(rbin.eff2.p, rbin.eff4.p, rbin.eff6.p, rbin.eff8.p, ncol = 2)
dev.off()

####################################
## Different reference years
####################################

# 1995: R1-95 - R8-95

r.mod1.tab95 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p + fnrb_ch95 + (1 | uniqid2), data = region.dat)
r.mod2.tab95 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch95 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod3.tab95 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod4.tab95 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod5.tab95 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p + fnrb_ch95 + (1 | uniqid2), data = region.dat)
r.mod6.tab95 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch95 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod7.tab95 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod8.tab95 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1995_p  + (ESSlrscale2 | uniqid2), data = region.dat)

screenreg(list(r.mod1.tab95, r.mod2.tab95, r.mod3.tab95, r.mod4.tab95, r.mod5.tab95, r.mod6.tab95, r.mod7.tab95, r.mod8.tab95), digits = 3)

rtab1 <- texreg(list(r.mod1.tab95, r.mod2.tab95, r.mod3.tab95, r.mod4.tab95), 
                file = "regionregs1_1995.tex", 
                stars = 0.05, 
                custom.model.names = c("R1-95", "R2-95", "R3-95", "R4-95"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1995)", "D Foreigners (1995-present)", "Left-Right x D Foreigners (1995-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results, using 1995 as base year. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs1_1995", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab2 <- texreg(list(r.mod5.tab95, r.mod6.tab95, r.mod7.tab95, r.mod8.tab95), 
                file = "regionregs2_1995.tex", 
                stars = 0.05, 
                custom.model.names = c("R5-95", "R6-95", "R7-95", "R8-95"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1995)", "D Foreigners (1995-present)", "Left-Right x D Foreigners (1995-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                center = TRUE, 
                caption = "Region-level results, using 1995 as base year. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs2_1995", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

# 2000: R1-00 - R8-00

r.mod1.tab00 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p + fnrb_ch00 + (1 | uniqid2), data = region.dat)
r.mod2.tab00 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch00 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod3.tab00 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod4.tab00 <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod5.tab00 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p + fnrb_ch00 + (1 | uniqid2), data = region.dat)
r.mod6.tab00 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch00 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod7.tab00 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod8.tab00 <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb2000_p  + (ESSlrscale2 | uniqid2), data = region.dat)

screenreg(list(r.mod1.tab00, r.mod2.tab00, r.mod3.tab00, r.mod4.tab00, r.mod5.tab00, r.mod6.tab00, r.mod7.tab00, r.mod8.tab00), digits = 3)

rtab1 <- texreg(list(r.mod1.tab00, r.mod2.tab00, r.mod3.tab00, r.mod4.tab00), 
                file = "/Users/johanneskarreth/Documents/Uni/1 - Papers/26 - Diversity by State/Data/Tables/regionregs1_2000.tex", 
                stars = 0.05, 
                custom.model.names = c("R1-00", "R2-00", "R3-00", "R4-00"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (2000)", "D Foreigners (2000-present)", "Left-Right x D Foreigners (2000-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                # reorder.coef = c(2, 13, 14, 15, 16, 17, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 18, 19, 1),
                center = TRUE, 
                caption = "Region-level results, using 2000 as base year. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs1_2000", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

rtab2 <- texreg(list(r.mod5.tab00, r.mod6.tab00, r.mod7.tab00, r.mod8.tab00), 
                file = "/Users/johanneskarreth/Documents/Uni/1 - Papers/26 - Diversity by State/Data/Tables/regionregs2_2000.tex", 
                stars = 0.05, 
                custom.model.names = c("R5-00", "R6-00", "R7-00", "R8-00"), 
                custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "Education", "Income (Decile)", "Unemployment", "% Foreigners (2000)", "D Foreigners (2000-present)", "Left-Right x D Foreigners (2000-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
                custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                digits = 3, 
                leading.zero = TRUE, 
                omit.coef = "(Intercept)", 
                # reorder.coef = c(2, 13, 14, 15, 16, 17, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 18, 19, 1),
                center = TRUE, 
                caption = "Region-level results, using 2000 as base year. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                caption.above = TRUE, 
                label = "tab:regionregs2_2000", 
                dcolumn = TRUE, 
                booktabs = TRUE, 
                use.packages = FALSE, 
                table = TRUE)

####################################
## Alternative immigration attitudes as outcome variables: R9 - R14
####################################

## Immigration bad (low) or good (high) for economy?
r.mod9.tab <- lmer(ESSimbgeco ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod10.tab <- lmer(ESSimbgeco ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod11.tab <- lmer(ESSimbgeco ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod12.tab <- lmer(ESSimbgeco ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)

## Immigration undermines (low) or enriches (high) country's cultural life
r.mod13.tab <- lmer(ESSimueclt ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod14.tab <- lmer(ESSimueclt ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod15.tab <- lmer(ESSimueclt ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod16.tab <- lmer(ESSimueclt ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)

## Immigrants make country worse (low) or better (high) place to live
r.mod17.tab <- lmer(ESSimwbcnt ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod18.tab <- lmer(ESSimwbcnt ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod19.tab <- lmer(ESSimwbcnt ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod20.tab <- lmer(ESSimwbcnt ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)

texreg(list(r.mod10.tab, r.mod12.tab, r.mod14.tab, r.mod16.tab, r.mod18.tab, r.mod20.tab), 
          file = "otherDVs.tex", 
          stars = 0.1, 
          custom.model.names = c("R9", "R10", "R11", "R12", "R13", "R14"), 
          custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Left/Right", "D Foreigners (1991-present)", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "Left-Right x D Foreigners (1991-present)", "% Non-Western Foreigners", "Left/Right x % Non-Western Foreigners"),
          custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
          digits = 3, 
          leading.zero = TRUE, 
          omit.coef = "(Intercept)", 
          reorder.coef = c(1, 2, 3, 5, 6, 7, 8, 4, 9, 10, 11),
          center = TRUE, 
          caption = "Region-level results, alternative immigration attitude items. Higher values on the outcome variables indicate higher acceptance.", 
          caption.above = TRUE, 
          label = "tab:otherDVs", 
          dcolumn = TRUE, 
          booktabs = TRUE, 
          use.packages = FALSE, 
          table = TRUE)

####################################
## Split samples L/C/R (region-level): R1-Left - R7-Right
####################################

r.mod1.left.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 < 5, ])
r.mod1.center.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 == 5, ])
r.mod1.right.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 > 5, ])

r.mod3.left.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 < 5, ])
r.mod3.center.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 == 5, ])
r.mod3.right.tab <- lmer(ESSimdfetn2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 > 5, ])

r.mod5.left.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 < 5, ])
r.mod5.center.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 == 5, ])
r.mod5.right.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 > 5, ])

r.mod7.left.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 < 5, ])
r.mod7.center.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 == 5, ])
r.mod7.right.tab <- lmer(ESSimpcntr2 ~ ESSfemale + ESSagegroup + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr + (1 | uniqid2), data = region.dat[region.dat$ESSlrscale2 > 5, ])

screenreg(list(r.mod1.left.tab, r.mod1.center.tab, r.mod1.right.tab, r.mod3.left.tab, r.mod3.center.tab, r.mod3.right.tab), digits = 3)

screenreg(list(r.mod5.left.tab, r.mod5.center.tab, r.mod5.right.tab, r.mod7.left.tab, r.mod7.center.tab, r.mod7.right.tab), digits = 3)

rtab1_lcr <- texreg(list(r.mod1.left.tab, r.mod1.center.tab, r.mod1.right.tab, r.mod3.left.tab, r.mod3.center.tab, r.mod3.right.tab), 
                    file = "regionregs1_lcr.tex", 
                    stars = 0.05, 
                    custom.model.names = c("R1-Left", "R1-Center", "R1-Right", "R3-Left", "R3-Center", "R3-Right"), 
                    custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "% Non-Western Foreigners"),
                    custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                    digits = 3, 
                    leading.zero = TRUE, 
                    omit.coef = "(Intercept)", 
                    center = TRUE, 
                    caption = "Region-level results, samples split in respondents on the political Left, Center, and Right. Outcome: ``Should your country allow many immigrants of different race/ethnic group from majority?'' Higher values on the outcome variables indicate higher acceptance.", 
                    caption.above = TRUE, 
                    label = "tab:regionregs1_lcr", 
                    dcolumn = TRUE, 
                    booktabs = TRUE, 
                    use.packages = FALSE, 
                    table = TRUE)

rtab2_lcr <- texreg(list(r.mod5.left.tab, r.mod5.center.tab, r.mod5.right.tab, r.mod7.left.tab, r.mod7.center.tab, r.mod7.right.tab), 
                    file = "regionregs2_lcr.tex", 
                    stars = 0.05, 
                    custom.model.names = c("R5-Left", "R5-Center", "R5-Right", "R7-Left", "R7-Center", "R7-Right"), 
                    custom.coef.names = c("Intercept", "Female", "Age (Decile)", "Education", "Income (Decile)", "Unemployment", "% Foreigners (1991)", "D Foreigners (1991-present)", "% Non-Western Foreigners"),
                    custom.note = "%stars. Cells: coefficients from hierarchical models with random intercepts for regions (standard errors in parentheses).", 
                    digits = 3, 
                    leading.zero = TRUE, 
                    omit.coef = "(Intercept)", 
                    center = TRUE, 
                    caption = "Region-level results, samples split in respondents on the political Left, Center, and Right. Outcome: ``Should your country allow many immigrants from poorer countries outside Europe?'' Higher values on the outcome variables indicate higher acceptance.", 
                    caption.above = TRUE, 
                    label = "tab:regionregs2_lcr", 
                    dcolumn = TRUE, 
                    booktabs = TRUE, 
                    use.packages = FALSE, 
                    table = TRUE)

## Plots

r.mod1.left.eff <- effect("fnrb_ch91", r.mod1.left.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnrb_ch91", layout = c(3, 1))
r.mod1.center.eff <- effect("fnrb_ch91", r.mod1.center.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)
r.mod1.right.eff <- effect("fnrb_ch91", r.mod1.right.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)

r.mod3.left.eff <- effect("nwepr", r.mod3.left.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnrb_ch91", layout = c(3, 1))
r.mod3.center.eff <- effect("nwepr", r.mod3.center.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)
r.mod3.right.eff <- effect("nwepr", r.mod3.right.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)

r.mod5.left.eff <- effect("fnrb_ch91", r.mod5.left.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnrb_ch91", layout = c(3, 1))
r.mod5.center.eff <- effect("fnrb_ch91", r.mod5.center.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)
r.mod5.right.eff <- effect("fnrb_ch91", r.mod5.right.tab, xlevels = list(fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)

r.mod7.left.eff <- effect("nwepr", r.mod7.left.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90); # plot(c.eff2, as.table = TRUE, x.var = "fnrb_ch91", layout = c(3, 1))
r.mod7.center.eff <- effect("nwepr", r.mod7.center.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)
r.mod7.right.eff <- effect("nwepr", r.mod7.right.tab, xlevels = list(nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)

r.mod1.left.eff.dat <- data.frame(r.mod1.left.eff$fit, r.mod1.left.eff$x, r.mod1.left.eff$lower, r.mod1.left.eff$upper, LR = "Left"); names(r.mod1.left.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod1.center.eff.dat <- data.frame(r.mod1.center.eff$fit, r.mod1.center.eff$x, r.mod1.center.eff$lower, r.mod1.center.eff$upper, LR = "Center"); names(r.mod1.center.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod1.right.eff.dat <- data.frame(r.mod1.right.eff$fit, r.mod1.right.eff$x, r.mod1.right.eff$lower, r.mod1.right.eff$upper, LR = "Right"); names(r.mod1.right.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod1.lcr.eff.dat <- rbind(r.mod1.left.eff.dat, r.mod1.center.eff.dat, r.mod1.right.eff.dat)

r.mod3.left.eff.dat <- data.frame(r.mod3.left.eff$fit, r.mod3.left.eff$x, r.mod3.left.eff$lower, r.mod3.left.eff$upper, LR = "Left"); names(r.mod3.left.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod3.center.eff.dat <- data.frame(r.mod3.center.eff$fit, r.mod3.center.eff$x, r.mod3.center.eff$lower, r.mod3.center.eff$upper, LR = "Center"); names(r.mod3.center.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod3.right.eff.dat <- data.frame(r.mod3.right.eff$fit, r.mod3.right.eff$x, r.mod3.right.eff$lower, r.mod3.right.eff$upper, LR = "Right"); names(r.mod3.right.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod3.lcr.eff.dat <- rbind(r.mod3.left.eff.dat, r.mod3.center.eff.dat, r.mod3.right.eff.dat)

r.mod5.left.eff.dat <- data.frame(r.mod5.left.eff$fit, r.mod5.left.eff$x, r.mod5.left.eff$lower, r.mod5.left.eff$upper, LR = "Left"); names(r.mod5.left.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod5.center.eff.dat <- data.frame(r.mod5.center.eff$fit, r.mod5.center.eff$x, r.mod5.center.eff$lower, r.mod5.center.eff$upper, LR = "Center"); names(r.mod5.center.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod5.right.eff.dat <- data.frame(r.mod5.right.eff$fit, r.mod5.right.eff$x, r.mod5.right.eff$lower, r.mod5.right.eff$upper, LR = "Right"); names(r.mod5.right.eff.dat) <- c("fit", "fnrb_ch91", "lower", "upper", "LR")
r.mod5.lcr.eff.dat <- rbind(r.mod5.left.eff.dat, r.mod5.center.eff.dat, r.mod5.right.eff.dat)

r.mod7.left.eff.dat <- data.frame(r.mod7.left.eff$fit, r.mod7.left.eff$x, r.mod7.left.eff$lower, r.mod7.left.eff$upper, LR = "Left"); names(r.mod7.left.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod7.center.eff.dat <- data.frame(r.mod7.center.eff$fit, r.mod7.center.eff$x, r.mod7.center.eff$lower, r.mod7.center.eff$upper, LR = "Center"); names(r.mod7.center.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod7.right.eff.dat <- data.frame(r.mod7.right.eff$fit, r.mod7.right.eff$x, r.mod7.right.eff$lower, r.mod7.right.eff$upper, LR = "Right"); names(r.mod7.right.eff.dat) <- c("fit", "nwepr", "lower", "upper", "LR")
r.mod7.lcr.eff.dat <- rbind(r.mod7.left.eff.dat, r.mod7.center.eff.dat, r.mod7.right.eff.dat)

rug.dat <- data.frame(fnrb_ch91 = region.dat$fnrb_ch91, nwepr = region.dat$nwepr)

r.mod1.lcr.eff.p <- ggplot(dat = r.mod1.lcr.eff.dat, aes(x = fnrb_ch91, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (R1-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

r.mod3.lcr.eff.p <- ggplot(dat = r.mod3.lcr.eff.dat, aes(x = nwepr, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("% Non−Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (R3-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

r.mod5.lcr.eff.p <- ggplot(dat = r.mod5.lcr.eff.dat, aes(x = fnrb_ch91, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (R5-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

r.mod7.lcr.eff.p <- ggplot(dat = r.mod7.lcr.eff.dat, aes(x = nwepr, y = fit)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper), alpha = 0.5) + geom_line() + facet_wrap(~ LR) + theme_bw() + guides (colour = FALSE, fill = FALSE) + xlab("% Non−Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (R7-Left/Center/Right)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") # +  scale_y_continuous(breaks = c(2.45 - 0.82, 2.45 - 0.5 * 0.82, 2.45, 2.45 + 0.5 * 0.82, 2.45 + 0.82), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD"))

pdf("int_region_lcr.pdf", width = 12, height = 10)
grid.arrange(r.mod1.lcr.eff.p, r.mod3.lcr.eff.p, r.mod5.lcr.eff.p, r.mod7.lcr.eff.p, ncol = 2)
dev.off()

####################################
## Figure 5: Conditional effects (region-level)
####################################

r.mod1 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod2 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod3 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod4 <- lmer(ESSimdfetn2.std ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod5 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + fnrb_ch91 + (1 | uniqid2), data = region.dat)
r.mod6 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup + ESSlrscale2*fnrb_ch91 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)
r.mod7 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup + ESSlrscale2 + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p + nwepr  + (1 | uniqid2), data = region.dat)
r.mod8 <- lmer(ESSimpcntr2.std ~ ESSfemale + ESSagegroup + ESSlrscale2*nwepr + ESSedulvla2 + ESShinc2 + unemp + fnrb1991_p  + (ESSlrscale2 | uniqid2), data = region.dat)

r.eff2 <- effect("ESSlrscale2*fnrb_ch91", r.mod2, xlevels = list(ESSlrscale2 = c(0, 5, 10), fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90); # plot(r.eff2, as.table = TRUE, x.var = "fnr_ch91", layout = c(3, 1))
r.eff4 <- effect("ESSlrscale2*nwepr", r.mod4, xlevels = list(ESSlrscale2 = c(0, 5, 10), nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)
r.eff6 <- effect("ESSlrscale2*fnrb_ch91", r.mod6, xlevels = list(ESSlrscale2 = c(0, 5, 10), fnrb_ch91 = sort(unique(region.dat$fnrb_ch91))), confidence.level = 0.90)
r.eff8 <- effect("ESSlrscale2*nwepr", r.mod8, xlevels = list(ESSlrscale2 = c(0, 5, 10), nwepr = sort(unique(region.dat$nwepr))), confidence.level = 0.90)

r.eff2.dat <- cbind(r.eff2$fit, r.eff2$x, r.eff2$lower, r.eff2$upper); names(r.eff2.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
r.eff2.dat$lr2 <- as.factor(r.eff2.dat$lr); levels(r.eff2.dat$lr2) <- c("Left", "Center", "Right")
r.eff4.dat <- cbind(r.eff4$fit, r.eff4$x, r.eff4$lower, r.eff4$upper); names(r.eff4.dat) <- c("fit", "lr", "nwepr", "lower", "upper")
r.eff4.dat$lr2 <- as.factor(r.eff4.dat$lr); levels(r.eff4.dat$lr2) <- c("Left", "Center", "Right")
r.eff6.dat <- cbind(r.eff6$fit, r.eff6$x, r.eff6$lower, r.eff6$upper); names(r.eff6.dat) <- c("fit", "lr", "fnrb_ch91", "lower", "upper")
r.eff6.dat$lr2 <- as.factor(r.eff6.dat$lr); levels(r.eff6.dat$lr2) <- c("Left", "Center", "Right")
r.eff8.dat <- cbind(r.eff8$fit, r.eff8$x, r.eff8$lower, r.eff8$upper); names(r.eff8.dat) <- c("fit", "lr", "nwepr", "lower", "upper")
r.eff8.dat$lr2 <- as.factor(r.eff8.dat$lr); levels(r.eff8.dat$lr2) <- c("Left", "Center", "Right")

rug.dat <- data.frame(fnrb_ch91 = region.dat$fnrb_ch91, nwepr = region.dat$nwepr)

r.eff2.p <- ggplot(dat = r.eff2.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (R2)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

r.eff4.p <- ggplot(dat = r.eff4.dat, aes(x = nwepr, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Ethnic acceptance (R4)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

r.eff6.p <- ggplot(dat = r.eff6.dat, aes(x = fnrb_ch91, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = fnrb_ch91, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("Change in % foreign (1991-present)") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (R6)") + geom_rug(data = rug.dat, aes(x = jitter(fnrb_ch91), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 2.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 2.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 2.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue")) # + theme(legend.position = c(0.1, 0.2))

r.eff8.p <- ggplot(dat = r.eff8.dat, aes(x = nwepr, y = fit, colour = lr2, group = lr2)) + geom_point(size = 0) + geom_ribbon(aes(x = nwepr, ymin = lower, ymax = upper, group = lr2, fill = lr2), alpha = 0.5) + geom_line() + theme_bw() + labs(fill = "Ideology") + guides (colour = FALSE, fill = FALSE) + xlab("% Non-Western foreigners") + ylab("Acceptance (low to high)") + ggtitle("Economic acceptance (R8)") + geom_rug(data = rug.dat, aes(x = jitter(nwepr), y = NULL, colour = NULL, group = NULL), sides = "b") + annotate("text", x = 7.5, y = 0.85, label = "Left", colour = "red") + annotate("text", x = 7.5, y = -0.2, label = "Center", colour = "gray") + annotate("text", x = 7.5, y = -0.85, label = "Right", colour = "blue") +  scale_y_continuous(breaks = c(-1, -0.5, 0, 0.5, 1), labels = c("-1 SD", "-0.5 SD", "Mean", "+0.5 SD", "+1 SD")) + scale_colour_manual(values=c("red", "gray", "blue")) + scale_fill_manual(values=c("red", "gray", "blue"))# + theme(legend.position = c(0.1, 0.2))

pdf("int_region.pdf", width = 8, height = 10)
grid.arrange(r.eff2.p, r.eff4.p, r.eff6.p, r.eff8.p, ncol = 2)
dev.off()

####################################
## Does ideology fully explain immigration attitudes?
####################################

summary(lm(ESSimdfetn2 ~ ESSlrscale2, data = region.dat))
cor(region.dat$ESSimdfetn2, region.dat$ESSlrscale2, use = "pairwise.complete.obs", method = "spearman")
dv1.tab <- table(region.dat$ESSimdfetn2, region.dat$ESSlrscale2)
polychor(dv1.tab, std.err = TRUE)
summary(lm(ESSimpcntr2 ~ ESSlrscale2, data = region.dat))
cor(region.dat$ESSimpcntr2, region.dat$ESSlrscale2, use = "pairwise.complete.obs", method = "spearman")
dv2.tab <- table(region.dat$ESSimpcntr2, region.dat$ESSlrscale2)
polychor(dv2.tab, std.err = TRUE)

p.dat <- melt(region.dat[, c("ESSimdfetn2", "ESSimpcntr2", "ESSlrscale2")], measure.vars = c("ESSimdfetn2", "ESSimpcntr2"))
p.dat <- p.dat[is.na(p.dat$value) == FALSE & is.na(p.dat$ESSlrscale2) == FALSE, ]
levels(p.dat$variable) <- c("Ethnic", "Economic")
p.dat <- mutate(group_by(p.dat, ESSlrscale2, value), obs.in.bin = n())
label.dat <- data.frame(ESSlrscale2 = c(3, 3), value = c(1.6, 1.4), label = c("tau~(Ethnic~acceptance) == 0.24", "tau~(Economic~acceptance) == 0.23"), variable = NA)
p1 <- ggplot(data = p.dat, aes(y = value, x = ESSlrscale2, colour = variable, fill = variable)) + geom_smooth() + geom_point(aes(alpha = 1/log(obs.in.bin + 1)), size = 1, position = position_jitter(w = 0.2, h = 0.2)) + guides(alpha = FALSE) + labs(fill = "Acceptance", colour = "Acceptance") + ylab("Acceptance (low to high)") + xlab("Ideology (left to right)") + theme_bw() + geom_text(data = label.dat, aes(x = ESSlrscale2, y = value, label = label), parse = TRUE, alpha = 1, colour = "black", size = 4) + scale_x_continuous(breaks = c(0, 5, 10))

png("ideology_acceptance.png", width = 8, height = 6, units = "in", res = 150)
p1
dev.off()

pdf("ideology_acceptance.pdf", width = 8, height = 6)
p1
dev.off()

####################################
## Summary statistics
####################################

library(pastecs)
library(knitr)

# Individual-level variables

summary.dat <- rbind(t(stat.desc(r.mod1.ryre@frame)), t(stat.desc(r.mod5.ryre@frame)))
summaryind.dat <- summary.dat[c("ESSimdfetn2", "ESSimpcntr2", "ESSfemale", "ESSagegroup", "ESSlrscale2", "ESSedulvla2", "ESShinc2"), c( "mean", "std.dev", "min", "max", "nbr.val")]
summaryind.tab <- round(summaryind.dat, digits = 2)
rownames(summaryind.tab) <- c("Acceptance (ethnic)", "Acceptance (economic)", "Female", "Age group", "Left/Right Self-Placement", "Edcuation", "Income (Decile)")
colnames(summaryind.tab) <- c("Mean/Proportion", "Std. dev.", "Min.", "Max.", "N (respondents)")
kable(summaryind.tab, format = "markdown", row.names = TRUE, caption = "Summary statistics of individual-level variables")

# Region-level variables

summaryreg1.dat <- summarise(group_by(r.mod1.ryre@frame, uniqid2yr), unemp = mean(unemp), fnrb1991_p = mean(fnrb1991_p), fnrb_ch91 = mean(fnrb_ch91))
summaryreg3.dat <- summarise(group_by(r.mod3.ryre@frame, uniqid2yr), nwepr = mean(nwepr))
summaryreg.dat <- data.frame(summaryreg1.dat, nwepr = summaryreg3.dat$nwepr)
summaryreg.dat <- t(stat.desc(summaryreg.dat))
summaryreg.dat <- summaryreg.dat[c("unemp", "fnrb1991_p", "fnrb_ch91", "nwepr"), c( "mean", "std.dev", "min", "max", "nbr.val")]
summaryreg.tab <- round(summaryreg.dat, digits = 2)
rownames(summaryreg.tab) <- c("Unemployment", "% Foreign (1991)", "∆ % Foreign (1991-present)", "% Non-Western")
colnames(summaryreg.tab) <- c("Mean/Proportion", "Std. dev.", "Min.", "Max.", "N (respondents)")
kable(summaryreg.tab, format = "markdown", row.names = TRUE, caption = "Summary statistics of region-level variables")